DESeq2 analyses steps
This section uses the counts data generated in Section 1 to differential expression analyses using DESeq2.
Prerequisites
Packages required for this analyses are loaded as follows:
setwd("~/github/amnion.vs.other_RNASeq")
# load the modules
library(sva)
library(tidyverse)
library(DESeq2)
library(vsn)
library(pheatmap)
library(ggrepel)
library(RColorBrewer)
library(reshape2)
require(biomaRt)
library(EnhancedVolcano)
library(TissueEnrich)
library(plotly)
library(DT)
library(cowplot)
library(biomaRt)Import datasets
The counts data and its associated metadata (coldata) are imported for analyses
counts = '~/github/amnion.vs.other_RNASeq/assets/counts-subset-v4.txt'
groupFile = '~/github/amnion.vs.other_RNASeq/assets/batch-subset-v4.txt'
coldata <- read.csv(groupFile, row.names=1, sep="\t", stringsAsFactors = TRUE)
cts <- as.matrix(read.csv(counts,sep="\t",row.names="gene.ids"))visualize the coldata
DT::datatable(coldata)order and inspect if all the coldata matches the headers of the counts data.
colnames(cts)
#> [1] "Naive_H9_hESCs_1" "Naive_H9_hESCs_2"
#> [3] "nTE_D1.Naive_H9_hESCs_1" "nTE_D1.Naive_H9_hESCs_2"
#> [5] "nTE_D2.Naive_H9_hESCs_1" "nTE_D2.Naive_H9_hESCs_2"
#> [7] "nTE_D3.Naive_H9_hESCs_1" "nTE_D3.Naive_H9_hESCs_2"
#> [9] "nCT_P3.Naive_H9_hESCs_1" "nCT_P3.Naive_H9_hESCs_2"
#> [11] "nCT_P10.Naive_H9_hESCs_1" "nCT_P10.Naive_H9_hESCs_2"
#> [13] "nCT_P15.Naive_H9_hESCs_1" "nCT_P15.Naive_H9_hESCs_2"
#> [15] "cR_nCT_P15.Naive_H9_hESCs_1" "cR_nCT_P15.Naive_H9_hESCs_2"
#> [17] "nCT_P15.409B2_iPSC_hESCs_1" "nCT_P15.409B2_iPSC_hESCs_2"
#> [19] "Placenta.derived_tbSCs_CT30_Ex1" "Placenta.derived_tbSCs_CT30_Ex2"
#> [21] "nST.Naive_H9_Ex1" "nST.Naive_H9_Ex2"
#> [23] "nEVT.Naive_H9_Ex1" "nEVT.Naive_H9_Ex2"
#> [25] "Primed_H9_hESCs_1" "Primed_H9_hESCs_2"
#> [27] "pBAP_D1.Primed_H9_hESCs_1" "pBAP_D1.Primed_H9_hESCs_2"
#> [29] "pBAP_D2.Primed_H9_hESCs_1" "pBAP_D2.Primed_H9_hESCs_2"
#> [31] "pBAP_D3.Primed_H9_hESCs_1" "pBAP_D3.Primed_H9_hESCs_2"
#> [33] "CytoTB_7_gestational_wks_1" "CytoTB_7_gestational_wks_2"
#> [35] "CytoTB_9_gestational_wks_1" "CytoTB_11_gestational_wks_1"
#> [37] "hESC_H1_STB_gt70um_D8_BAP_1" "hESC_H1_STB_gt70um_D8_BAP_2"
#> [39] "hESC_H1_STB_gt70um_D8_BAP_3" "hESC_H1_STB_40.70um_D8_BAP_1"
#> [41] "hESC_H1_STB_40.70um_D8_BAP_2" "hESC_H1_STB_40.70um_D8_BAP_3"
#> [43] "hESC_H1_STB_lt40um_D8_BAP_1" "hESC_H1_STB_lt40um_D8_BAP_2"
#> [45] "hESC_H1_STB_lt40um_D8_BAP_3" "hESC_H1_D8_MEF.CM.and.FGF2_1"
#> [47] "hESC_H1_D8_MEF.CM.and.FGF2_2" "hESC_H1_D8_MEF.CM.and.FGF2_3"
all(rownames(coldata) %in% colnames(cts))
#> [1] TRUE
cts <- cts[, rownames(coldata)]Batch correction
We will use the combat seq (SVA package) to run batch correction based on the bioproject (dataset origin)
cov1 <- as.factor(coldata$BioProject)
adjusted_counts <- ComBat_seq(cts, batch=cov1, group=NULL)
#> Found 2 batches
#> Using null model in ComBat-seq.
#> Adjusting for 0 covariate(s) or covariate level(s)
#> Estimating dispersions
#> Fitting the GLM model
#> Shrinkage off - using GLM estimates for parameters
#> Adjusting the data
all(rownames(coldata) %in% colnames(cts))
#> [1] TRUE
cts <- cts[, rownames(coldata)]run DESeq2
The batch corrected read counts are then used for running DESeq2 analyses
dds <- DESeqDataSetFromMatrix(countData = adjusted_counts,
colData = coldata,
design = ~ condition)
vsd <- vst(dds, blind=FALSE)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
dds
#> class: DESeqDataSet
#> dim: 16496 48
#> metadata(1): version
#> assays(4): counts mu H cooks
#> rownames(16496): ENSG00000000003.15 ENSG00000000005.6 ...
#> ENSG00000288695.1 ENSG00000288698.1
#> rowData names(106): baseMean baseVar ... deviance maxCooks
#> colnames(48): cR_nCT_P15.Naive_H9_hESCs_1 cR_nCT_P15.Naive_H9_hESCs_2
#> ... hESC_H1_STB_lt40um_D8_BAP_2 hESC_H1_STB_lt40um_D8_BAP_3
#> colData names(3): BioProject condition sizeFactorVarious contrasts are setup as follows (a total of 9 combinations)
res.nCTvsBAP <- results(dds, contrast=c("condition",
"nCT_P10.Naive_H9_hESCs",
"pBAP_D3.Primed_H9_hESCs"))
res.nTEvsSTB <- results(dds, contrast=c("condition",
"nTE_D3.Naive_H9_hESCs",
"hESC_H1_STB_gt70um_D8_BAP"))
res.nCTvsSTB <- results(dds, contrast=c("condition",
"nCT_P10.Naive_H9_hESCs",
"hESC_H1_STB_gt70um_D8_BAP"))
res.nTEvsBAP <- results(dds, contrast=c("condition",
"nTE_D3.Naive_H9_hESCs",
"pBAP_D3.Primed_H9_hESCs"))
res.nTEvsL40 <- results(dds, contrast=c("condition",
"nTE_D3.Naive_H9_hESCs",
"hESC_H1_STB_lt40um_D8_BAP"))
res.nCTvsL40 <- results(dds, contrast=c("condition",
"nCT_P10.Naive_H9_hESCs",
"hESC_H1_STB_lt40um_D8_BAP"))
res.BAPvsL40 <- results(dds, contrast=c("condition",
"pBAP_D3.Primed_H9_hESCs",
"hESC_H1_STB_lt40um_D8_BAP"))
res.BAPvsSTB <- results(dds, contrast=c("condition",
"pBAP_D3.Primed_H9_hESCs",
"hESC_H1_STB_gt70um_D8_BAP"))
res.STBvsL40 <- results(dds, contrast=c("condition",
"hESC_H1_STB_gt70um_D8_BAP",
"hESC_H1_STB_lt40um_D8_BAP"))function to save results
processDE <- function(restable, fileName) {
restable <- restable[order(restable$padj), ]
outTable <- merge(as.data.frame(restable),
as.data.frame(counts(dds, normalized=TRUE)),
by="row.names", sort=FALSE)
names(outTable)[1] <- "Gene"
newName=paste0("assets/DESeq2results-", fileName, "_fc.tsv")
write_delim(outTable, file=newName, delim = "\t")
upName <- paste0(fileName, ".up")
upName <-outTable %>% filter(log2FoldChange >= 1) %>%
filter(padj <= 0.05) %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(Gene)
downName <- paste0(fileName, ".donw")
downName <-outTable %>% filter(log2FoldChange <= 1) %>%
filter(padj <= 0.05) %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(Gene)
}The results are saved as tsv files
processDE(res.nCTvsBAP, "nCTvsBAP")
processDE(res.nCTvsSTB, "nCTvsSTB")
processDE(res.nTEvsSTB, "nTEvsSTB")
processDE(res.nTEvsBAP, "nTEvsBAP")
processDE(res.nTEvsL40, "nTEvsL40")
processDE(res.nCTvsL40, "nCTvsL40")
processDE(res.BAPvsL40, "BAPvsL40")
processDE(res.STBvsL40, "STBvsL40")
processDE(res.BAPvsSTB, "BAPvsSTB")Preparing GeneLists for PCE
We will create a function to save the DE gene lists (up and down) and run PCE on them we will need the annotations to convert the gene-id-version to just gene-id (and add other metadata later)
ensembl=useMart("ENSEMBL_MART_ENSEMBL")
listDatasets(ensembl) %>%
filter(str_detect(description, "Human"))
#> dataset description version
#> 1 hsapiens_gene_ensembl Human genes (GRCh38.p13) GRCh38.p13
ensembl = useDataset("hsapiens_gene_ensembl", mart=ensembl)
listFilters(ensembl) %>%
filter(str_detect(name, "ensembl"))
#> name
#> 1 ensembl_gene_id
#> 2 ensembl_gene_id_version
#> 3 ensembl_transcript_id
#> 4 ensembl_transcript_id_version
#> 5 ensembl_peptide_id
#> 6 ensembl_peptide_id_version
#> 7 ensembl_exon_id
#> description
#> 1 Gene stable ID(s) [e.g. ENSG00000000003]
#> 2 Gene stable ID(s) with version [e.g. ENSG00000000003.15]
#> 3 Transcript stable ID(s) [e.g. ENST00000000233]
#> 4 Transcript stable ID(s) with version [e.g. ENST00000000233.10]
#> 5 Protein stable ID(s) [e.g. ENSP00000000233]
#> 6 Protein stable ID(s) with version [e.g. ENSP00000000233.5]
#> 7 Exon ID(s) [e.g. ENSE00000000003]
filterType <- "ensembl_gene_id_version"
filterValues <- rownames(cts)
listAttributes(ensembl) %>%
head(20)
#> name description
#> 1 ensembl_gene_id Gene stable ID
#> 2 ensembl_gene_id_version Gene stable ID version
#> 3 ensembl_transcript_id Transcript stable ID
#> 4 ensembl_transcript_id_version Transcript stable ID version
#> 5 ensembl_peptide_id Protein stable ID
#> 6 ensembl_peptide_id_version Protein stable ID version
#> 7 ensembl_exon_id Exon stable ID
#> 8 description Gene description
#> 9 chromosome_name Chromosome/scaffold name
#> 10 start_position Gene start (bp)
#> 11 end_position Gene end (bp)
#> 12 strand Strand
#> 13 band Karyotype band
#> 14 transcript_start Transcript start (bp)
#> 15 transcript_end Transcript end (bp)
#> 16 transcription_start_site Transcription start site (TSS)
#> 17 transcript_length Transcript length (including UTRs and CDS)
#> 18 transcript_tsl Transcript support level (TSL)
#> 19 transcript_gencode_basic GENCODE basic annotation
#> 20 transcript_appris APPRIS annotation
#> page
#> 1 feature_page
#> 2 feature_page
#> 3 feature_page
#> 4 feature_page
#> 5 feature_page
#> 6 feature_page
#> 7 feature_page
#> 8 feature_page
#> 9 feature_page
#> 10 feature_page
#> 11 feature_page
#> 12 feature_page
#> 13 feature_page
#> 14 feature_page
#> 15 feature_page
#> 16 feature_page
#> 17 feature_page
#> 18 feature_page
#> 19 feature_page
#> 20 feature_page
attributeNames <- c('ensembl_gene_id_version',
'ensembl_gene_id',
'external_gene_name')
annot <- getBM(attributes=attributeNames,
filters = filterType,
values = filterValues,
mart = ensembl)
# remove duplicates, if any
isDup <- duplicated(annot$ensembl_gene_id)
dup <- annot$ensembl_gene_id[isDup]
annot <- annot[!annot$ensembl_gene_id%in%dup,]function to create up regulated gene list
pceUP <- function(restable) {
restable <- restable[order(restable$padj), ]
outTable <- merge(as.data.frame(restable),
as.data.frame(counts(dds, normalized=TRUE)),
by="row.names", sort=FALSE)
names(outTable)[1] <- "Gene"
upName <-outTable %>% filter(log2FoldChange >= 1) %>%
filter(padj <= 0.05) %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(Gene)
upNew <- annot[annot$ensembl_gene_id_version %in% upName$Gene,]
upList <- upNew$ensembl_gene_id
return(upList)
}function to create down regulated gene list
pceDOWN <- function(restable) {
restable <- restable[order(restable$padj), ]
outTable <- merge(as.data.frame(restable),
as.data.frame(counts(dds, normalized=TRUE)),
by="row.names", sort=FALSE)
names(outTable)[1] <- "Gene"
downName <-outTable %>% filter(log2FoldChange <= -1) %>%
filter(padj <= 0.05) %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(Gene)
downNew <- annot[annot$ensembl_gene_id_version %in% downName$Gene,]
downList <- downNew$ensembl_gene_id
return(downList)
}Run the functions on the DESeq2 results object
nCTvsBAP.up <- pceUP(res.nCTvsBAP)
nCTvsSTB.up <- pceUP(res.nCTvsSTB)
nTEvsSTB.up <- pceUP(res.nTEvsSTB)
nTEvsBAP.up <- pceUP(res.nTEvsBAP)
nTEvsL40.up <- pceUP(res.nTEvsL40)
nCTvsL40.up <- pceUP(res.nCTvsL40)
BAPvsL40.up <- pceUP(res.BAPvsL40)
STBvsL40.up <- pceUP(res.STBvsL40)
BAPvsSTB.up <- pceUP(res.BAPvsSTB)
nCTvsBAP.down <- pceDOWN(res.nCTvsBAP)
nCTvsSTB.down <- pceDOWN(res.nCTvsSTB)
nTEvsSTB.down <- pceDOWN(res.nTEvsSTB)
nTEvsBAP.down <- pceDOWN(res.nTEvsBAP)
nTEvsL40.down <- pceDOWN(res.nTEvsL40)
nCTvsL40.down <- pceDOWN(res.nCTvsL40)
BAPvsL40.down <- pceDOWN(res.BAPvsL40)
STBvsL40.down <- pceDOWN(res.STBvsL40)
BAPvsSTB.down <- pceDOWN(res.BAPvsSTB)Run PCE
The above gene lists are used for running PCE. The function used for running PCE is below
# load the PCE data
l <- load(file = "~/TutejaLab/PlacentaEnrich/combine-test-expression1.Rdata")
humanGeneMapping <- dataset$GRCH38$humanGeneMapping
d <- dataset$PlacentaDeciduaBloodData
data <- d$expressionData
cellDetails <- d$cellDetails
# create a run PCE function
runpce <- function(inputgenelist, title, shade) {
inputGenes<-toupper(inputgenelist)
expressionData<-data[intersect(row.names(data),humanGeneMapping$Gene),]
se<-SummarizedExperiment(assays = SimpleList(as.matrix(expressionData)),
rowData = row.names(expressionData),
colData = colnames(expressionData))
cellSpecificGenesExp<-teGeneRetrieval(se,expressedGeneThreshold = 1)
print(length(inputGenes))
gs<-GeneSet(geneIds=toupper(inputGenes))
output2<-teEnrichmentCustom(gs,cellSpecificGenesExp)
enrichmentOutput<-setNames(data.frame(assay(output2[[1]]),
row.names = rowData(output2[[1]])[,1]),
colData(output2[[1]])[,1])
row.names(cellDetails)<-cellDetails$RName
enrichmentOutput$Tissue<- cellDetails[row.names(enrichmentOutput),"CellName"]
ggplot(data = enrichmentOutput, mapping = aes(x = reorder (Tissue, -Log10PValue), Log10PValue)) +
geom_bar(stat = "identity", color = shade, fill = shade) + theme_classic(base_size = 10) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 10),
plot.title = element_text(size = 12),
plot.margin = unit(c(1,1,1,2), "cm")) +
labs(x = "",
y = "-log10 p-value (adj.)",
title = title) +
scale_y_continuous(expand = expansion(mult = c(0, .1)))
}The PCE is run on each of the gene list as follows (up and down pair are displayed together)
p <- runpce(nCTvsBAP.up , "overexpressed in H9_nCT_P10_Io", "#F16746")
#> [1] 1381
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nCTvsBAP.down , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 2821
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nCT_P10_Io vs. BAP.d3
p <- runpce(nCTvsSTB.up , "overexpressed in H9_nCT_P10_Io)", "#F16746")
#> [1] 563
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nCTvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 1766
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nCT_P10_Io vs. H1_BAP_D8_>70_Yabe
p <- runpce(nTEvsSTB.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 1069
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nTEvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 2070
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nTE_D3_Io vs. H1_BAP_D8_>70_Yabe
p <- runpce(nTEvsBAP.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 886
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nTEvsBAP.down , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 2160
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nTE_D3_Io vs. pBAP_D3_Io
p <- runpce(nTEvsL40.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 1183
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nTEvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 2260
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nTE_D3_Io vs. H1_BAP_D8_<40_Yabe
p <- runpce(nCTvsL40.up , "overexpressed in H9_nCT_P10_Io", "#F16746")
#> [1] 650
#> No background list provided. Using all the
#> genes as background.
q <- runpce(nCTvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 2095
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H9_nCT_P10_Io vs. H1_BAP_D8_<40_Yabe
p <- runpce(BAPvsL40.up , "up regulated in pBAP_D3_Io", "#0773B2")
#> [1] 1544
#> No background list provided. Using all the
#> genes as background.
q <- runpce(BAPvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 1351
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe
p <- runpce(STBvsL40.up , "overexpressed in H1_BAP_D8_>70_Yabe", "#AFBF38")
#> [1] 1266
#> No background list provided. Using all the
#> genes as background.
q <- runpce(STBvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#0773B2")
#> [1] 1157
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: H1_BAP_D8_>70_Yabe vs. H1_BAP_D8_<40_Yabe
p <- runpce(BAPvsSTB.up , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 1732
#> No background list provided. Using all the
#> genes as background.
q <- runpce(BAPvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 1424
#> No background list provided. Using all the
#> genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plotFigX: pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe
Volcano plots
We will setup a funciton for drawing volcano plots and then run them on the DESeq2 results. The funciton is as shown below
mart <- read.csv("~/TutejaLab/amnion_for_ms_20210715/de-analyses/mart-genes.tsv",
sep="\t", stringsAsFactors = TRUE,
header = TRUE)
myVolPlot <- function(restable, first, second, shade1, shade2) {
restable <- restable[order(restable$padj), ]
outTable <- tibble::rownames_to_column(as.data.frame(restable), "Gene")
outTable <- merge(outTable, mart, by.x=c("Gene"),
by.y = c("ensembl_gene_id_version"),
sort=FALSE)
outTable <- outTable %>%
mutate_all(na_if,"") %>%
mutate_all(na_if," ") %>%
mutate(gene_symbol = coalesce(gene_symbol,Gene))
outTable$diffexpressed <- "other.genes"
outTable$diffexpressed[outTable$log2FoldChange >= 1 &
outTable$padj <= 0.05] <- paste0("Higher expression in ", first)
outTable$diffexpressed[outTable$log2FoldChange <= -1 &
outTable$padj <= 0.05] <- paste0("Higher expression in ", second)
outTable$delabel <- ""
outTable$delabel[outTable$log2FoldChange >= 1
& outTable$padj <= 0.05
& !is.na(outTable$padj)] <- outTable$gene_symbol[outTable$log2FoldChange >= 1
& outTable$padj <= 0.05
& !is.na(outTable$padj)]
outTable$delabel[outTable$log2FoldChange <= -1
& outTable$padj <= 0.05
& !is.na(outTable$padj)] <- outTable$gene_symbol[outTable$log2FoldChange <= -1
& outTable$padj <= 0.05
& !is.na(outTable$padj)]
ggplot(outTable, aes(x=log2FoldChange, y=-log10(padj), col=diffexpressed, label=delabel)) +
geom_point(alpha = 0.5) +
theme_classic() +
scale_color_manual(name = "Expression", values=c(shade1, shade2, "#4d4d4d")) +
# ggtitle(paste0(first, " vs. ", second)) +
xlab("log2 fold change") +
ylab("-log10 pvalue (adjusted)") +
theme(legend.text.align = 0,
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(5, 5, 5, 5))
}Running Volcano plots for each comparison is shown below
ggplotly(myVolPlot(res.nCTvsBAP, "H9_nCT_P10_Io", "H9_pBAP_D3_Io", "#F16746", "#0571b0" ))FigX: H9_nCT_P10_Io vs. H9_pBAP_D3_Io
ggplotly(myVolPlot(res.nCTvsSTB, "H9_nCT_P10_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#F16746" ))FigX: H9_nCT_P10_Io vs. H1_BAP_D8_>70_Yabe
ggplotly(myVolPlot(res.nTEvsSTB, "H9_nTE_D3_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#C79D2E" ))FigX: H9_nTE_D3_Io vs. H1_BAP_D8_>70_Yabe
ggplotly(myVolPlot(res.nTEvsBAP, "H9_nTE_D3_Io", "H9_pBAP_D3_Io", "#C79D2E", "#0571b0" ))FigX: H9_nTE_D3_Io vs. H9_pBAP_D3_Io
ggplotly(myVolPlot(res.nTEvsL40, "H9_nTE_D3_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#C79D2E"))FigX: H9_nTE_D3_Io vs. H1_BAP_D8_<40_Yabe
ggplotly(myVolPlot(res.nCTvsL40, "H9_nCT_P10_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#F16746" ))FigX: H9_nCT_P10_Io vs. H1_BAP_D8_<40_Yabe
ggplotly(myVolPlot(res.BAPvsL40, "H9_pBAP_D3_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#0571b0" ))FigX: H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe
ggplotly(myVolPlot(res.STBvsL40, "H1_BAP_D8_>70_Yabe", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#5C823A" ))FigX: H1_BAP_D8_>70_Yabe vs. H1_BAP_D8_<40_Yabe
ggplotly(myVolPlot(res.BAPvsSTB, "H9_pBAP_D3_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#0571b0" ))FigX: H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe